Generation and analysis of SymEnergy single storage model results

Companion notebook #1 of 4 of the paper

M.C.Soini et al., On the market displacement of incumbent grid-connected electricity storage by more efficient storage.

This notebook runs the single-storage SymEnergy model and compare the model's results with manually constructed results as shown in figure 4 of the paper.

Instructions

  • Install dependencies as defined in the README

Imports

In [1]:
import sys
import pandas as pd
import numpy as np
import sympy as sp
import itertools
from tqdm import tqdm
from collections import namedtuple
import time
from bokeh.io import show, output_notebook
output_notebook(verbose=False)
sp.init_printing(pretty_print=False)
Loading BokehJS ...

SymEnergy imports

In [2]:
from symenergy import Model
from symenergy import Evaluator
from symenergy import logger
from symenergy.evaluator import plotting
from symenergy import cache_params
from symenergy import multiproc_params
cache_params['path'] = './symenergy_cache/single_storage'
logger.setLevel('WARNING')

Four time slot model, single POWER-constrained storage

In [3]:
# Additional assumption on the behavior of the system:
# g (dispatchable power plant) power production non-zero 
# -> all positivity constraints of dispatchable 
#    power production are non-binding ("False") 
constraint_filt = ('act_lb_g_pos_p_0nig == False and '
                   'act_lb_g_pos_p_1mor == False and '
                   'act_lb_g_pos_p_2day == False and '
                   'act_lb_g_pos_p_3eve == False'
                  )

m_pwr = Model(curtailment=False, nworkers=10, constraint_filt=constraint_filt)

m_pwr.add_slot(name='0nig', load=10000, weight=6)
m_pwr.add_slot(name='1mor', load=10000, weight=6)
m_pwr.add_slot(name='2day', load=10000, weight=6)
m_pwr.add_slot(name='3eve', load=10000, weight=6)

m_pwr.add_plant(name='g', vc0=0, vc1=1)
m_pwr.add_storage(name='phs', eff=0.99, capacity=1, energy_cost=1e-12)

m_pwr.generate_solve()
> 15:45:25 - WARNING - symenergy.auxiliary.io - ****************************************************************
> 15:45:25 - WARNING - symenergy.auxiliary.io - ****************************************************************
> 15:45:25 - WARNING - symenergy.auxiliary.io - Loading from cache file .../single_storage/m59E4716525F6.pickle.
> 15:45:25 - WARNING - symenergy.auxiliary.io - Please delete this file to re-solve model: Model.cache.delete()
> 15:45:25 - WARNING - symenergy.auxiliary.io - ****************************************************************
> 15:45:25 - WARNING - symenergy.auxiliary.io - ****************************************************************

Evaluation

Evaluate the model for selected numerical parameter values and plot the results.

In [4]:
logger.setLevel('WARNING')
multiproc_params['nworkers'] = 'default'

x_vals = {
          m_pwr.slots['0nig'].l: np.linspace(3000, 10000, 7),
          m_pwr.slots['1mor'].l: np.linspace(3000, 10000, 7),
          m_pwr.slots['2day'].l: np.linspace(3000, 10000, 7),
          m_pwr.slots['3eve'].l: np.linspace(9000, 10000, 2),
          m_pwr.comps['phs'].eff: np.linspace(0.6, 0.999, 3),
          m_pwr.comps['phs'].C: np.linspace(0, 1000, 3),
         }
ev = Evaluator(m_pwr, x_vals, drop_non_optimum=True)
ev.get_evaluated_lambdas_parallel()
multiproc_params['nworkers'] = 1
ev.expand_to_x_vals_parallel()
> 15:45:34 - WARNING - symenergy.auxiliary.io - **********************************************************************
> 15:45:34 - WARNING - symenergy.auxiliary.io - **********************************************************************
> 15:45:34 - WARNING - symenergy.auxiliary.io - Loading from cache file .../single_storage/l4FDD524A2300.pickle.
> 15:45:34 - WARNING - symenergy.auxiliary.io - Please delete this file to re-evaluate: Evaluator.cache_lambd.delete()
> 15:45:34 - WARNING - symenergy.auxiliary.io - **********************************************************************
> 15:45:34 - WARNING - symenergy.auxiliary.io - **********************************************************************
> 15:45:34 - WARNING - symenergy.auxiliary.io - *********************************************************************
> 15:45:34 - WARNING - symenergy.auxiliary.io - *********************************************************************
> 15:45:34 - WARNING - symenergy.auxiliary.io - Loading from cache file .../single_storage/e7B92BC306C05.pickle.
> 15:45:34 - WARNING - symenergy.auxiliary.io - Please delete this file to re-evaluate: Evaluator.cache_eval.delete()
> 15:45:34 - WARNING - symenergy.auxiliary.io - *********************************************************************
> 15:45:34 - WARNING - symenergy.auxiliary.io - *********************************************************************

Interactive plot of model results as a function of the PHS capacity

In [5]:
balplot = plotting.BalancePlot(ev, ind_axx='C_phs_none', ind_pltx='slot', ind_plty=None, plot_width=300)
show(balplot._get_layout())
BokehDeprecationWarning: 'WidgetBox' is deprecated and will be removed in Bokeh 3.0, use 'bokeh.models.Column' instead

Four time slot model, single ENERGY-constrained storage

In [6]:
logger.setLevel('WARNING')
multiproc_params['nworkers'] = 'default'

# g (dispatchable power plant) power production non-zero -> all g positivity constraint non-binding ("False") 
constraint_filt = ('act_lb_g_pos_p_0nig == False and '
                   'act_lb_g_pos_p_1mor == False and '
                   'act_lb_g_pos_p_2day == False and '
                   'act_lb_g_pos_p_3eve == False')

m_erg = Model(curtailment=False, nworkers=10, constraint_filt=constraint_filt)

m_erg.add_slot(name='0nig', load=10000, weight=6)
m_erg.add_slot(name='1mor', load=10000, weight=6)
m_erg.add_slot(name='2day', load=10000, weight=6)
m_erg.add_slot(name='3eve', load=10000, weight=6)

m_erg.add_plant(name='g', vc0=0, vc1=1)
m_erg.add_storage(name='phs', eff=0.99, energy_capacity=1, energy_cost=1e-12, charging_to_energy_factor='eta')

m_erg.generate_solve()
> 15:45:35 - WARNING - symenergy.auxiliary.io - ****************************************************************
> 15:45:35 - WARNING - symenergy.auxiliary.io - ****************************************************************
> 15:45:35 - WARNING - symenergy.auxiliary.io - Loading from cache file .../single_storage/mF1B02BF3B06F.pickle.
> 15:45:35 - WARNING - symenergy.auxiliary.io - Please delete this file to re-solve model: Model.cache.delete()
> 15:45:35 - WARNING - symenergy.auxiliary.io - ****************************************************************
> 15:45:35 - WARNING - symenergy.auxiliary.io - ****************************************************************

Evaluation

Evaluate the model for certain numerical parameter values and plot the results.

In [7]:
logger.setLevel('WARNING')
multiproc_params['nworkers'] = 'default'

x_vals = {
          m_erg.slots['0nig'].l: np.linspace(3000, 10000, 7),
          m_erg.slots['1mor'].l: np.linspace(3000, 10000, 7),
          m_erg.slots['2day'].l: np.linspace(3000, 10000, 7),
          m_erg.slots['3eve'].l: np.linspace(9000, 10000, 2),
          m_erg.comps['phs'].eff: [0.9],
          m_erg.comps['phs'].E: np.linspace(0, 6000, 4),
         }
ev = Evaluator(m_erg, x_vals, drop_non_optimum=False)
ev.get_evaluated_lambdas_parallel()
multiproc_params['nworkers'] = 1
ev.expand_to_x_vals_parallel()
> 15:45:38 - WARNING - symenergy.auxiliary.io - **********************************************************************
> 15:45:38 - WARNING - symenergy.auxiliary.io - **********************************************************************
> 15:45:38 - WARNING - symenergy.auxiliary.io - Loading from cache file .../single_storage/lFBA5F1A37B67.pickle.
> 15:45:38 - WARNING - symenergy.auxiliary.io - Please delete this file to re-evaluate: Evaluator.cache_lambd.delete()
> 15:45:38 - WARNING - symenergy.auxiliary.io - **********************************************************************
> 15:45:38 - WARNING - symenergy.auxiliary.io - **********************************************************************
> 15:45:38 - WARNING - symenergy.auxiliary.io - *********************************************************************
> 15:45:38 - WARNING - symenergy.auxiliary.io - *********************************************************************
> 15:45:38 - WARNING - symenergy.auxiliary.io - Loading from cache file .../single_storage/e48CCAE8CBE68.pickle.
> 15:45:38 - WARNING - symenergy.auxiliary.io - Please delete this file to re-evaluate: Evaluator.cache_eval.delete()
> 15:45:38 - WARNING - symenergy.auxiliary.io - *********************************************************************
> 15:45:38 - WARNING - symenergy.auxiliary.io - *********************************************************************
In [8]:
balplot = plotting.BalancePlot(ev, ind_axx='E_phs_none', ind_pltx='slot', ind_plty=None, plot_width=300)
show(balplot())
BokehDeprecationWarning: 'WidgetBox' is deprecated and will be removed in Bokeh 3.0, use 'bokeh.models.Column' instead

Comparison to constructed solutions

Verification of the blueprint in the upper part of figure 4 by comparison to solutions obtained from the Lagrange approach.

Definition of auxiliary functions

The following cells define various auxiliary functions for the manual construction of the results.

In [9]:
tolist = lambda lst: '(%s)' % ', '.join(lst)

def get_operation(x, model):
    ''' Get strings encoding operational patterns. '''
    
    batphs = 'phs'
    slots = model.slots
    is_E_capacity_constrained = 'E_phs_none' in model.comps['phs'].parameters('name')
    is_C_capacity_constrained = 'C_phs_none' in model.comps['phs'].parameters('name')
    
    list_c = ([False if x[f'act_lb_{batphs}_pos_pchg_{slot}'] else 'c' for slot in slots])
    list_d = ([False if x[f'act_lb_{batphs}_pos_pdch_{slot}'] else 'd' for slot in slots])

    if is_C_capacity_constrained:
        list_c = [c.upper() if x[f'act_lb_{batphs}_pchg_cap_C_{slot}'] and c else c for slot, c in zip(slots, list_c)]
        list_d = [d.upper() if x[f'act_lb_{batphs}_pdch_cap_C_{slot}'] and d else d for slot, d in zip(slots, list_d)]
    
    list_cd = ['-' if not c and not d else (c if c else d) for c, d in zip(list_c, list_d)]

    list_eC = (['E' if is_E_capacity_constrained and x[f'act_lb_{batphs}_e_cap_E_{slot}'] else False for slot in slots])
    list_e0 = (['0' if x[f'act_lb_{batphs}_pos_e_{slot}'] else False for slot in slots])

    list_e = ['-' if not c and not d else (c if c else d) for c, d in zip(list_eC, list_e0)]
    
    idling = list_cd == ['-'] * 4
    name_batphs = batphs
    return pd.Series({f'cd_pattern_{name_batphs}':  ' ' + tolist(list_e) + ' '+ tolist(list_cd),
                      f'is_idling_{name_batphs}': idling})
In [10]:
def get_blocks(df, model):
    '''
    Get blocks of time slots separated by zero and/or capacity-constrained energy
    '''
    
    slots = model.slots
    is_E_capacity_constrained = 'E_phs_none' in model.comps['phs'].parameters('name')

    if is_E_capacity_constrained:
        list_cap_E = df[[f'act_lb_phs_e_cap_E_{slot}' for slot in slots]].astype(int).tolist()
    else:
        list_cap_E = [0, 0, 0, 0]
        
    list_pos_E = df[[f'act_lb_phs_pos_e_{slot}' for slot in slots]].astype(int).tolist()

    if sum(list_cap_E) == 0 and sum(list_pos_E) == 1:
        return [0] * 4

    df_blocks = pd.DataFrame({'act_cap_E': list_cap_E, 'act_pos_E': list_pos_E, 'slot_id': range(4)})

    df_blocks = pd.concat([df_blocks.iloc[[-1], :],
                           df_blocks])

    df_blocks['diff_cap'] = df_blocks.act_cap_E.diff()
    df_blocks['diff_pos'] = df_blocks.act_pos_E.diff()

    if df_blocks.act_cap_E.any() and -1 in df_blocks['diff_cap'].values:
        shift = np.where(df_blocks['diff_cap'].values == -1)[0][0] - 1 
    elif df_blocks.act_pos_E.sum() > 0 and -1 in df_blocks['diff_pos'].values:
        shift = np.where(df_blocks['diff_pos'].values == -1)[0][0] - 1
    else: shift = 0

    
    df_blocks = df_blocks.iloc[1:]
    df_blocks['slot_id_shift'] = np.roll(df_blocks.slot_id, int(shift))
    df_blocks = df_blocks.sort_values('slot_id_shift')

    df_blocks = pd.concat([df_blocks.iloc[[-1], :],
                           df_blocks])

    df_blocks['diff_cap'] = df_blocks.act_cap_E.diff()
    df_blocks['diff_pos'] = df_blocks.act_pos_E.diff()

    block = [0] *4
    iblock = 0
    is_discharging = True
    for i, (cap, pos, last_cap, last_pos) in enumerate(list(zip(df_blocks.act_cap_E, df_blocks.act_pos_E, 
                 df_blocks.act_cap_E.shift(1), df_blocks.act_pos_E.shift(1)))[1:]):

        if last_pos == 1:
            is_discharging = False
        if last_cap == 1:
            is_discharging = True

        start_after_cap = cap == 0 and last_cap == 1
        start_after_pos = pos == 0 and last_pos == 1
        stays_full = cap == 1 and last_cap == 1
        stays_empty = pos == 1 and last_pos == 1

        if start_after_cap or start_after_pos or stays_full or stays_empty:
            iblock += 1

        block[i] = iblock * (10 if is_discharging else 1)

    return np.roll(block, shift)

Get manual solution components

Various functions to get the components as shown in figure 4 of the paper ("unconstrained operation drivers", "power constraint terms", etc)

In [11]:
def get_single_N(this_charging, other_charging, this_slot, other_slot):
    other_w = getattr(w, other_slot)
    this_l = getattr(l, this_slot)
    other_l = getattr(l, other_slot)

    return {(True, True):  other_w * eff**2 * (other_l - this_l),
            (True, False):  other_w * (eff * other_l - this_l),
            (False, True):  other_w * eff * (other_l - this_l * eff),
            (False, False):  other_w * (other_l - this_l),
            }[(this_charging, other_charging)]

def get_single_D(other_charging, other_slot):
    other_w = getattr(w, other_slot)

    return {(True): other_w * eff**2,
            (False): other_w
           }[other_charging]

def get_single_N_C(this_charging, other_charging, other_slot):
    other_w = getattr(w, other_slot)

    return {(True, True):  - other_w * C * eff**2,
            (True, False): + other_w * C * eff,
            (False, True): - other_w * C * eff,
            (False, False):  + other_w * C,
            }[(this_charging, other_charging)]

def get_N_E(this_charging, is_charging_block):

    return {(True, True):  E * eff,
            (True, False): E,
            (False, True): - E * eff,
            (False, False): - E,
            }[(is_charging_block, this_charging)]
    
def get_sign(this_charging):
    
    return {True: +1, False: -1}[this_charging]

def get_N_unconstr(is_charging, slots, slct_idx, same_block_indices, is_active, is_power_cap_constr):
    return sum(get_single_N(this_charging=is_charging[slct_idx], 
                            other_charging=is_charging[other_idx],
                            this_slot=slots[slct_idx],
                            other_slot=slots[other_idx])
               for other_idx in same_block_indices
               if is_active[other_idx] and not is_power_cap_constr[other_idx]) 


def get_D(is_charging, slots, same_block_indices, is_active, is_power_cap_constr):
    return sum(get_single_D(other_charging=is_charging[other_idx],
                            other_slot=slots[other_idx])
               for other_idx in same_block_indices
               if is_active[other_idx] and not is_power_cap_constr[other_idx])

def get_N_C(is_charging, slots, slct_idx, same_block_indices, is_active, is_power_cap_constr):
    return sum(get_single_N_C(this_charging=is_charging[slct_idx],
                              other_charging=is_charging[other_idx],
                              other_slot=slots[other_idx])
            for other_idx in same_block_indices
            if is_active[other_idx] and is_power_cap_constr[other_idx])

The compare_row function

For each row (constraint combination) of the solution table, construct the result manually and compare with the model's result.

In [12]:
def compare_row(row, model):

    is_E_capacity_constrained = 'E_phs_none' in model.comps['phs'].parameters('name')
    is_C_capacity_constrained = 'C_phs_none' in model.comps['phs'].parameters('name')
    slots = [slot[1:] for slot in model.slots]
    nslots = list(model.slots)
    
    if row.is_idling_phs:
        return
        
    is_charging = (1 - row[[r for r in row.index if 'pos_pchg' in r]].astype(int)).values
    is_discharging = (1 - row[[r for r in row.index if 'pos_pdch' in r]].astype(int)).values
    is_cap_E_constrained = (row[[r for r in row.index if 'cap_E' in r]].astype(int).values
                            if is_E_capacity_constrained
                            else np.zeros(4))
    is_power_cap_constr = (np.array([row[f'act_lb_phs_pchg_cap_C_{slot}'] or row[f'act_lb_phs_pdch_cap_C_{slot}'] for slot in model.slots])
                           if is_C_capacity_constrained
                           else np.zeros(4))
    is_zero_E = row[[r for r in row.index if 'pos_e' in r]].astype(int).values

    is_active = np.array([(a or b) for a, b in zip(is_charging, is_discharging)])
    block_indices = np.array(row.blocks)
    (slct_indices,) = np.where(is_active)
    any_binding_E = row[[r for r in row.index if 'cap_E' in r]].sum() > 0

    # loop over slots with active storage power
    for slct_idx in slct_indices:
        (same_block_indices,) = np.where(np.array(block_indices) == block_indices[slct_idx])
        is_charging_block = block_indices[same_block_indices][0] < 10

        is_incomplete_block = False

        # last index *before* the block and last index of the block
        dict_prev_last_ind = {
         (0, 1): (3, 1), (1, 2): (0, 2), (2, 3): (1, 3), (3, 0): (2, 0),
         (0, 1, 2): (3, 2), (1, 2, 3): (0, 3), (2, 3, 0): (1, 0), (3, 0, 1): (2, 1)}
        dict_prev_last_ind = {tuple(set(k)): v for k, v in dict_prev_last_ind.items()}
        ind_before_block, last_ind = (dict_prev_last_ind[tuple(set(same_block_indices))] 
                                      if tuple(set(same_block_indices)) in dict_prev_last_ind 
                                      else (None, None))

        # partial discharging/charging between two zero E/two capacity constrained slots
        if ind_before_block is not None:
            is_incomplete_block = (len(same_block_indices) > 1
                                   and ((is_cap_E_constrained[ind_before_block] and is_cap_E_constrained[last_ind])
                                        or (is_zero_E[ind_before_block] and is_zero_E[last_ind])))

        # get constructed solution
        if not is_power_cap_constr[slct_idx]:
            N = get_N_unconstr(is_charging, slots, slct_idx, same_block_indices, is_active, is_power_cap_constr)
            N += get_N_C(is_charging, slots, slct_idx, same_block_indices, is_active, is_power_cap_constr)
            if any_binding_E and not is_incomplete_block:
                N += get_N_E(this_charging=is_charging[slct_idx], is_charging_block=is_charging_block)
            D = get_D(is_charging, slots, same_block_indices, is_active, is_power_cap_constr)
            sgn = get_sign(this_charging=is_charging[slct_idx])
            result_constr = sgn * N / D
        else:
            result_constr = C

        # get solution from Lagrange approach
        var_name = f'phs_p{"chg" if is_charging[slct_idx] else "dch"}_{nslots[slct_idx]}'
        result_expect = model.get_results_dict(idx=row.idx, slct_var_mlt=[var_name], substitute={'ec_phs_none': 0})[var_name]
        dict_subs_powers = {2.0: 2, eff**1.0: eff, 3.0: 3, 4.0: 4, 5.0: 5, 6.0: 6, 7.0: 7, 8.0: 8}
        result_expect = sp.simplify(sp.simplify(result_expect.subs(dict_subs_powers)).subs(dict_subs_powers))

        # compare solutions; 
        if sp.simplify(result_expect - result_constr) != 0:
            print('idx: ', row.idx, 'slot:',  slct_idx)
            print('EXPECT', result_expect)
            print('CONSTR', result_constr)
            print('slct_idx', slct_idx, slots)
            print('is_power_cap_constr', is_power_cap_constr)
            print('cd_pattern_phs', row.cd_pattern_phs)

            raise ValueError('Solutions don\'t match')

Auxiliary function to obtain the sympy symbols

In [13]:
def get_sympy_symbols(model):
    '''Retrieve the sympy symbols from the model instance
    required to construct the model solutions.'''
    
    slots = np.array(['nig', 'mor', 'day', 'eve'])
    nslots = list(map(''.join, zip(map(str, range(4)), slots)))
    L = namedtuple('L', slots)
    W = namedtuple('W', slots)

    # extract symbols from model parameters
    (vc1_g_none, 
     l_0nig, w_0nig, l_1mor, w_1mor, 
     l_2day, w_2day, l_3eve, w_3eve,
     ec_phs_none, eff_phs_none, 
     E_phs_none, vre_scale_none
    ) = model.parameters('symb')

    # aggregate for convenience
    l = L(l_0nig, l_1mor, l_2day, l_3eve)
    w = W(w_0nig, w_1mor, w_2day, w_3eve)
    eff = eff_phs_none
    E = E_phs_none
    
    return l, w, eff, E

Verify all solutions of both models

Applied the compare_row function to each row of the result DataFrame. If no errors are thrown, all constructed results are identical to the ones obtained from the model.get_results_dict method.

In [14]:
for model, pwrerg in [(m_erg, 'energy'), (m_pwr, 'power ')]:

    l, w, eff, E = get_sympy_symbols(model)
    C = E

    df_comb = pd.concat([model.df_comb, model.df_comb.apply(get_operation, model=model, axis=1)], axis=1)

    tqdm.pandas(desc=f'Get blocks for model with {pwrerg} constraints      ')
    df_comb['blocks'] = df_comb.progress_apply(get_blocks, model=model, axis=1)

    tqdm.pandas(desc=f'Compare solutions of model with {pwrerg} constraints')
    df_comb.progress_apply(compare_row, model=model, axis=1)
/home/martinsoini/.conda/envs/symenergy_final_single/lib/python3.8/site-packages/tqdm/std.py:666: FutureWarning: The Panel class is removed from pandas. Accessing it from the top-level namespace will also be removed in the next version
  from pandas import Panel
Get blocks for model with energy constraints      : 100%|██████████| 121/121 [00:00<00:00, 240.99it/s]
Compare solutions of model with energy constraints: 100%|██████████| 121/121 [01:10<00:00,  1.71it/s]
Get blocks for model with power  constraints      : 100%|██████████| 463/463 [00:00<00:00, 670.97it/s]
Compare solutions of model with power  constraints: 100%|██████████| 463/463 [03:28<00:00,  2.22it/s]